femR is an R package that solves PDEs on R relying on
the Finite Element Method (FEM). Let \(\Omega
\subset \mathbb{R}^d\) be the domain of our interest,
femR package solves advection-reaction-diffusion problem of
the following form:
\[
\begin{cases}
-\mu \Delta f + \boldsymbol{b} \cdot \nabla f + c f = u \qquad & in
\ \Omega \\
BC \ conditions \ (Dirichlet \ or \ Neumann) \qquad & on \
\partial \Omega,
\end{cases}
\] where \(\mu \in \mathbb{R}\),
\(\boldsymbol{b} \in \mathbb{R}^d\),
\(c \in \mathbb{R}\) are the diffusion
coefficient, the transport coefficient and the reaction coefficient,
respectively. This vignette illustrates how to use femR to
solve a partial differential equation. In particular, the following
section explains step by step the solution of a simple Poisson
problem.
Let \(\Omega\) be the unit square \([0,1]^2\). Consider the following problem defined over \(\Omega\):
\[ \begin{cases} -\Delta f = u &in \ \Omega \\ \quad f = 0 &on \ \partial \Omega \end{cases} \]
where \(u(x,y) = 8\pi^2 \sin( 2\pi x)\sin(2\pi y)\) is the forcing term and \(\partial \Omega\) is the boundary of \(\Omega\) where we have prescribed homogeneous Dirichelet boundary condition. The exact solution of the previous problem is \(u_{ex}(x,y) = \sin(2\pi x) \sin(2 \pi y)\). The following figure shows the exact solution of the problem, on the left hand side, and the triangulation of the unit square domain, on the right hand side.
The following window wraps all the steps needed to solve the problem.
## 1. Defining domain
data("unit_square", package="femR")
mesh <- Mesh(unit_square)
## 2. Defining the solution of the PDE
Vh <- FunctionSpace(mesh, fe_order=2)
f <- f <- Function(Vh)
## 3. Defining the differential operator
L <- -laplace(f)
## 4. Defining the forcing term and the Dirichlet boundary conditions as standard R functions
# forcing term
u <- function(points){
return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2]))
}
# Dirichlet BC
dirichletBC <- function(points){
return(matrix(0,nrow=nrow(points), ncol=1))
}
## 5. Building the PDE object
pde <- Pde(L, u, dirichletBC)
## 6. computing the discrete solution
pde$solve()
## 7. Plots
plot(f)fig_h <- contour(f)
subplot(fig_exact %>%hide_colorbar(), fig_h, margin = 0.05) %>%layout(annotations = list(
list(x=0.155, y=1.05, text = "exact solution", showarrow = F,
xref='paper', yref='paper'),
list(x=0.875, y=1.05, text = "discrete solution", showarrow = F,
xref='paper', yref='paper')))The following sections explain more in detail the previous chunk of code.
First, we create a proper domain object exploiting the function
create_mesh which takes a R named list as parameter. Such
named list should contain geometrical information of the domain:
nodes, a 2-columns matrix containing the coordinates
of the nodes of the mesh.
elements, a 3-columns matrix that for each row \(i\) contains the indexes of the mesh nodes
that are vertices of the \(i\)-th
element.
boundary, a 1-column matrix that indicates the nodes
on the boundary of the domain.
femR contains the unit_square named list
that can be used to perform some tests.
data("unit_square", package="femR")
unit_square <- Mesh(unit_square)
plot(unit_square)Moreover, femR package provides an utility that reads
files storing such geometrical information and provided by third-party
software such as FreeFem++. Indeed,
read_mesh(filename) reads a .mesh file from
FreeFem++ and returns a named list that can be passed to
create_mesh function. The following chunk, that is not
executed, shows how read_mesh works.
filename <- "path/to/domain.mesh" #domain.mesh stored by the FreeFem++ utility savemesh()
domain <- read_mesh(filename)
mesh <- Mesh(domain)First, we define initialize a FunctionSpace object
passing the mesh tha we previously build and the finite element order.
Note that, the package provide first order and second order finite
elements. First order finite element is the default parameter. Then, we
define a Function, solution of the PDEs, belonging to the
function space Vh.
## solution of the PDE
Vh <- FunctionSpace(unit_square, fe_order=2)
f <- Function(Vh)Then, we define the differential operator in its strong formulation.
L <- -laplace(f) # L <- -div(I*grad(f))
# In general:
# L <- (-mu)*laplace(f) + dot(b,grad(f)) + c*fNote that, according to the well-known identity \(div(\nabla f) = \Delta f\), we could have
written L<--div(I*grad(f)), where I is the
2-dimensional identity matrix. Indeed, femR package let you
define every term of a generic second-order linear differential operator
\(Lf= -div(K\nabla f) + \mathbf{b}\cdot \nabla
f + a f\) as follows:
Diffusion: -div(K*grad(f)) or
(-mu)*laplacian(f), where either K or
mu represents the diffusion matrix in an anisotropic
problem or the diffusion coefficient in a isotropic problem, implements
the first term of the differential operator according to the problem at
hand.
Advection: dot(b,grad(f)), where b is a
two-dimensional vector, implements the second term of differential
operator.
Reaction: c*f, where c is a scalar,
implements the last term of the differential operator.
Although, the interface let the user define a generic second-order
linear differential operator: -div(K*grad(f)) + ..., at
this stage the package has been tested only in the isotropic case, i.e
K=I, where I is the two-dimensional identy
matrix.
We define the forcing term of the differential problem and the
Dirichlet boundary condition as a simple R functions. Such
functions will be provided as parameters to the object which
encapsulates the concept of PDEs.
## forcing term
u <- function(points){
return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2]))
}
## Dirichlet BC
dirichletBC <- function(points){
return(matrix(0,nrow=nrow(points), ncol=1))
}Finally, we build a pde object passing the differential
operator L, as first parameter, the forcing term
u, as second parameter, the Dirichlet boundary condition
dirichletBC, as third parameter:
## create pde
pde <- Pde(L, u, dirichletBC)We can compute the discrete solution of the Poisson problem calling
the solve method:
## solve problem
pde$solve()Since, we are dealing with a simple differential problem and we know
the analytic form of the exact solution, we can compute the \(L^2\) norm of the error exploiting on the
mass matrix provided by the pde object:
exact_solution <- function(points){
return( sin(2.*pi*points[,1])*sin(2.*pi*points[,2]))
}
u_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates()))
error.L2 <- sqrt(sum(pde$get_mass()%*%(u_ex-pde$solution())^2))
cat("L2 error = ", error.L2, "\n")
#> L2 error = 0.0004136347The R package plotly is exploited to plot the computed
solution.
## plot solution
plot(f)# contour
contour(f)Note that, the base::plot function and
graphics::contourfunction have both been overloaded and
return a plotly object that can be modified as one wishes.
plot(f, colorscale="Jet") %>% layout(scene=list(aspectmode="cube"))